#include "header.h"

static int N;
static FILE *fh_out;
static double *k, *v, g;
static fftw_complex *a[4], *RS;

void init_prep(params_ptr in) {
  N = in->n;
  v = fftw_malloc(sizeof(double)*N);
  k = fftw_malloc(sizeof(double)*N);
  for (int j = 0; j < in->n; j ++) {
    v[j] = 2*pi*(j+0.5)/N;
    k[j] = j;
    if (j > N/2) k[j] = j - N;
  }
  for (int j = 0; j < 4; j++) a[j] = fftw_malloc(N*sizeof(fftw_complex));
  RS = fftw_malloc(N*sizeof(fftw_complex));
  g = in->g;
}

void write_out(fftw_complex* in, char* str) {
  fh_out = fopen(str, "w");
  fprintf(fh_out, "#1.u 2.k 3.re R 4. im R\n\n");
  for (int j = 0; j < N; j++) {
    fprintf(fh_out, "%e\t%e\t%e\t%e\n", v[j], k[j], creal(in[j]), cimag(in[j]) );
  }
  fclose(fh_out);
}


void write_out_surface(work_ptr in, char* str) {
  fh_out = fopen(str, "w");
  fprintf(fh_out, "#1.u 2.x 3.y\n\n");
  r_transform(in->rho, RS, in->alpha);
  for (int j = 0; j < N; j++) {
    fprintf(fh_out, "%e\t%e\t%e\n", v[j], creal(v[j] + 1I*(in->alpha)+(in->rho[j]) + 1.I*RS[j])/(2*pi),
                                          cimag(v[j] + 1I*(in->alpha)+(in->rho[j]) + 1.I*RS[j])/(2*pi) );
  }
  fclose(fh_out);
}
